################################################################################
#
# Moon.py - A Python module that ...
#
# Copyright (C) 2011 Adrian Price-Whelan
# 
################################################################################

"""
todo: something to do
"""

__author__ = 'Adrian Price-Whelan <adrn@nyu.edu>'

# Standard library dependencies (e.g. sys, os)
import math
import datetime

# Internal imports
import convert
from Sun import Sun

class Moon(object):
	"""
	"""
	ra = None
	dec = None
	datetimeObj = None
	
	def __init__(self, datetimeObj=datetime.datetime.now()):
		self.datetimeObj = datetimeObj
		self._set_position()
	
	def _set_position(self):
		""" ~ a few arcminutes accuracy
		"""
		datetimeObj = self.datetimeObj
		
		l0 = 318.351648 # mean longitude
		P0 = 36.340410 # mean longitude of perigee
		N0 = 318.510107 # mean longitude of node
		ii = 5.145396 # inclination
		ee = 0.054900 # eccentricity
		aa = 384401 # km, semi-major axis or moon's orbit
		theta0 = 0.5181 # degrees, semiangular size at distance a
		pi0 = 0.9507 # parallax at distance a
		
		sun = Sun(datetimeObj)
		
		jdJan0 = convert.datetime2jd(datetime.datetime(datetimeObj.year,1,1,hour=0, minute=0, second=0))
		jd = convert.datetime2jd(datetimeObj)
		
		d = jd - jdJan0
		D = (datetimeObj.year - 1990) * 365.0 + (datetimeObj.year - 1992) / 4 + d + 2
		
		l = (13.1763966*D + l0) % 360.0
		C = l - sun.longitude

		moonMeanAnomaly = (l - 0.1114041*D - P0) % 360.0
		N = (N0 - 0.0529539*D) % 360.0
		
		Ev = 1.2739*math.sin(math.radians(2*C - moonMeanAnomaly))
		Ae = 0.1858*math.sin(math.radians(sun.meanAnomaly))
		A3 = 0.37 * math.sin(math.radians(sun.meanAnomaly))
				
		corrected_moonMeanAnomaly = moonMeanAnomaly + Ev - Ae - A3
		
		Ec = 6.2886 * math.sin(math.radians(corrected_moonMeanAnomaly))
		A4 = 0.214 * math.sin(math.radians(2.0 * corrected_moonMeanAnomaly))
		
		lprime = l + Ev + Ec - Ae + A4
		V = 0.6583 * math.sin(math.radians(2.0*(lprime - sun.longitude)))
		lprimeprime = lprime + V
		
		Nprime = N - 0.16*math.sin(math.radians(sun.meanAnomaly))
		y = math.sin(math.radians(lprimeprime - Nprime))*math.cos(math.radians(ii))
		x = math.cos(math.radians(lprimeprime - Nprime))
		
		arcTan = math.degrees(math.atan(y/x))
		
		if y > 0 and x > 0:
			arcTan = arcTan % 90.0
		elif y > 0 and x < 0:
			arcTan = (arcTan % 90.0) + 90.0
		elif y < 0 and x < 0:
			arcTan = (arcTan % 90.0) + 180.0
		elif y < 0 and x > 0:
			arcTan = (arcTan % 90.0) + 270.0
		
		moonLongitude = arcTan + Nprime
		moonBeta = math.degrees(math.asin(math.sin(math.radians(lprimeprime - Nprime)) * math.sin(math.radians(ii))))
		
		ra, dec = convert.eclipticLatLon2RADec(moonLongitude, moonBeta)
		
		self.ra = ra
		self.dec = dec
	
	def illumination(self, datetimeObj=datetime.datetime.now()):
		"""
		"""
		fraction = 0.0
		return fraction
	
	def rise(self, datetimeObj=datetime.datetime.now()):
		"""
		"""
		return datetimeObj
		
	def set(self, datetimeObj=datetime.datetime.now()):
		"""
		"""
		return datetimeObj
	

def main():
	moon1 = Moon(datetime.datetime(1979, 2, 26, 16, 00, 50))
	moon2 = Moon(datetime.datetime(2011, 3, 8, 3, 00, 00))
	print convert.dec2sex(moon2.ra/15.0), convert.dec2sex(moon2.dec)
	

if __name__ == '__main__':
	main()